function bin=iflt_simple(img, nbin)
if nargin < 2
     P = 8;
end
% if nargin < 3
%     R = 1;
% end
if nargin < 4
    nbin = 64;
end
if size(img,3) ~= 1
    img = rgb2gray(img);
end
img = double(img);
[M, N] = size(img);
% Ix = zeros(M,N,P);
% Iy = zeros(M,N,P);
np = ceil(P/2);
bin = zeros(nbin,4);
% maximum range for mean = -1 .. 1
% maximum range for std = 0..(np/(np-1))^0.5
rng_m = 1.0;
rng_s = (np/(np-1)) ^ 0.5;
%rng_s = (P/(P-1)) ^ 0.5;
u = 0.5*(nbin-1);
pm = u / rng_m; %scale factor
ps = u * 2 / rng_s; %scale factor

tic;
for i = 2:M-1
    for j = 2:N-1
        %Extract I
        a = img(i,j) - img(i-1:i+1,j-1:j+1);
        I = [a(1:4) a(6:9)];
        if ~any(I)
            f = round([u u u u])+1;
        else
            %Normalisation
            I_norm = I / ((sum(I.^2)^0.5)+eps);
            I_norm_1 = I_norm(1:2:end);
            I_norm_2 = I_norm(2:2:end);
            cA1 = (I_norm_1 + I_norm_2) / 2^0.5;
            cD1 = (I_norm_1 - I_norm_2) / 2^0.5;
            %f1->rng1 f2->rng2
            f1 = round([mean(cA1) mean(cD1)]*pm + u); % values in [0..nbin-1]
            f2 = round([std(cA1) std(cD1)]*ps); % values in [0..nbin-1]
            f = [f1 f2]+1;
        end
        s = sub2ind([nbin, 4], f, 1:4);
        bin(s) = bin(s) + 1;
    end
end
toc;
b=bin(:);
end
